Applying Sparse Matrix Solvers to a Glacial Ice Sheet Model

More recent models have attempted to use better physics to describe ice velocities in regions where velocities vary significantly over short distances. These models use a 3-dimensional, rectangular FEM grid. While the 2-dimensional model is used for an entire ice sheet, the 3-dimensional model is used over a limited area of interest with a smaller distance between grid points. Output of the 2-dimensional model is used to establish boundary conditions for the 3-dimensional model. When solving for 3D velocities, FEM generates systems of banded linear equations with 81 non-zero coefficients per equation. Unlike the 2D model, these equations are not diagonally dominant. For this reason, the decision was made to investigate direct methods based on Gaussian elimination for solving these systems of equations. Abstract Two software packages for solving sparse systems of linear equations, SuperLU and UMFPACK, have been integrated with the University of Maine Ice Sheet Model for predicting the formation and disappearance of glacial ice sheets. Using a library of basic linear algebra subprograms, BLAS, tuned for the underlying hardware, these packages perform significantly better than our banded Gaussian elimination routine. Also, they are able to solve non-banded systems that can be produced by the ice sheet model, but are impractical to solve with straightforward Gaussian elimination. A modified compressed column data structure is presented for interfacing the ice sheet model with the two software packages. Test results are presented that indicate careful consideration must be given to the column ordering methods used by the packages when solving specific problems. An additional complication in the 3D model is internal pressure within the ice sheet. In one version of the model the internal pressure is eliminated, resulting in a system of linear equations that is purely banded. In a second version of the 3D model, pressure is explicitly calculated and results in a banded system of equations with lower and right borders. Figure 1 depicts the non-zero entries in these systems of equations.